#!/bin/bash -e
# based on methy_anno_170122.sh, has no vcf in input
function info() {
echo Usage: `basename $0` 'methy_result field_chr,field_pos,field_beta,depth'
exit 2
}

while getopts  ":p:f:" opt; do
	case  $opt  in
		p) out_prefix=$OPTARG;;
		f) suffix=$OPTARG;;
		*) info;;
	esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi

. $var

tab2vcf="tab2vcf.sh -p$out_prefix"

methyl_id=methyl_id
beta_id=beta
depth_id=depth
bed450=$data_path/intervals/1/450.clean.$genome_assembly.bed.gz
vcf_header=$data_path/methy/vcf_header
beta=$out_prefix.beta.txt
# b10k=$data_path/methy/beta/b10k_0.txt
b3k=$data_path/methy/beta/b3k.txt

# Total size (bp): 80454470
# interval=$data/bed/capture.bed
site_total=$((5624218/2))


cat $1|cut -f$2|sed 1d > $beta
bgzip -c $beta > $beta.gz
tabix -f -b2 -e2 $beta.gz

$tab2vcf $beta 1,2

total=`grep -v '^#' $out_prefix.convert.vcf |wc -l`
echo -e "calling_total\t$total" > $out_prefix.stat.txt


#-------------------------filter-------------------------
filter_vcf_bed.sh -p$out_prefix.capture $out_prefix.convert.vcf $data_path/bed/capture.bed

vcf=$out_prefix.capture.filter.vcf
echo -e "site_total\t$site_total" >> $out_prefix.stat.txt
num_capture=`grep -v '^#' $vcf |wc -l` 
echo -e "site_capture\t$num_capture" >> $out_prefix.stat.txt
awk 'BEGIN{print "percent_capture\t" 100 * '"$num_capture"' / '"$site_total"'}' >> $out_prefix.stat.txt

#-----------------------------annotate beta value & depth-----------------
cat $vcf |sed /^##vcfProcessLog/d > $vcf.tmp && mv $vcf $vcf.tmp1 && mv $vcf.tmp $vcf && rm $vcf.tmp1
bgzip -c $vcf > $vcf.gz
tabix -pvcf $vcf.gz


echo beta annotate
bcftools annotate -a$beta.gz -h$vcf_header -o$out_prefix.beta.depth.vcf -Ov -c CHROM,POS,$beta_id,$depth_id $vcf.gz

bgzip -c $out_prefix.beta.depth.vcf > $out_prefix.beta.depth.vcf.gz
tabix -pvcf $out_prefix.beta.depth.vcf.gz

#-------------------------------filter and annotate to 450k---------------
echo methylation id annotate result region $bed450
bcftools annotate -R$bed450 -a$bed450 -h$vcf_header -o$out_prefix.cg450.vcf -Ov -c CHROM,FROM,TO,$methyl_id $out_prefix.beta.depth.vcf.gz
num_id=`grep -v '^#' $out_prefix.cg450.vcf |wc -l`
echo -e "num_450k\t$num_id" >> $out_prefix.stat.txt



echo extract
$java_run/snpsift extractField $out_prefix.cg450.vcf $methyl_id $beta_id|sed 1d > $out_prefix.cg450.beta.txt


# echo -e '\tbeta\ncancer\tna' |cat - $out_prefix.cg450.beta.txt > $out_prefix.cg450.beta.pd.txt



echo rf pre
rf_pre1.sh $out_prefix.cg450.beta.txt 10000

beta_stat.py $out_prefix.b450k.sub.test.txt $out_prefix.cg450.beta.txt
mv beta_mean $out_prefix.cg450.beta_mean.txt

rf_methy.sh $out_prefix.cg450.beta_mean.txt $out_prefix.b450k.sub.test.txt
rf_methy.sh $out_prefix.cg450.beta.txt $out_prefix.b450k.sub.test.txt

rf_methy.sh $out_prefix.cg450.beta_mean.txt $out_prefix.b450k.sub.test.sub.txt
rf_methy.sh $out_prefix.cg450.beta.txt $out_prefix.b450k.sub.test.sub.txt


# echo -e '\tbeta\ncancer\tna' |cat - beta_mean > $out_prefix.beta_mean.pd.txt

# rf_methy_011017.py $b10k $out_prefix.cg450.beta.pd.txt 200 4
# rf_methy_011017.py $b3k $out_prefix.cg450.beta.pd.txt 200 4
# rf_methy_011017.py $b3k $out_prefix.beta_mean.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.txt $out_prefix.cg450.beta.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.can_merge.txt $out_prefix.cg450.beta.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.can_sub.merge.txt $out_prefix.cg450.beta.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.txt $out_prefix.beta_mean.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.can_merge.txt $out_prefix.cg450.beta_mean.pd.txt 200 4

# rf_methy_011017.py $out_prefix.b450k.sub.test.can_sub.merge.txt $out_prefix.cg450.beta_mean.pd.txt 200 4

bk.sh -tproject/methylation/test/`date +%Y-%m` predict_score.txt

. $cmd_done